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Abstract 

The nature of freezing and melting transitions for a system of hard disks in 
a spatially periodic external potential is studied using extensive Monte Carlo 
simulations. Detailed finite size scaling analysis of various thermodynamic 
quantities like the order parameter, its cumulants etc. are used to map the 
phase diagram of the system for various values of the density and the ampli- 
tude of the external potential. We find clear indication of a re -entrant liquid 
phase over a significant region of the parameter space. Our simulations there- 
fore show that the system of hard disks behaves in a fashion similar to charge 
stabilized colloids which are known to undergo an initial freezing, followed 
by a re -melting transition as the amplitude of the imposed, modulating field 
produced by crossed laser beams is steadily increased. Detailed analysis of 
our data shows several features consistent with a recent dislocation unbinding 
theory of laser induced melting. 
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I. INTRODUCTION 



The liquid -solid transition in systems of particles under the influence of external mod- 
ulating potentials has recently attracted a fair amount of attention from experiments |l] {7| . 
theory |||| and computer simulations [p~0| p~3[1 . This is partly due to the fact that well con- 
trolled, clean experiments can be performed using colloidal particles \TA\ confined between 
glass plates (producing essentially a two -dimensional system) and subjected to a spatially 
periodic electromagnetic field generated by interfering two, or more, crossed laser beams. 
One of the more surprising results of these studies, where a commensurate, one dimensional, 
modulating potential is imposed, is the fact that there exist regions in the phase diagram 
over which one observes re -entrant freezing/melting behaviour. As a function of the 
laser field intensity the system first freezes from a modulated liquid to a two dimensional 
triangular solid - a further increase of the intensity confines the particles strongly within the 
troughs of the external potential, making the system quasi -one -dimensional which increases 
fluctuations and leads to re -melting. 

Our present understanding of this curious phenomenon has come from early mean -field 
density functional || and more recent dislocation unbinding || calculations. The mean field 
theories neglect fluctuations and therefore cannot explain re -entrant behaviour. The order 
of the transition is predicted to be first order for small laser field intensities, though for 
certain combinations of external potentials (which includes the specific geometry studied 
in the experiments and in this paper) the transition may become second order after going 
through a tricritical point. In general, though mean field theories are applicable in any 
dimension, the results are expected to be accurate only for higher dimensions and long 
ranged potentials. The validity of the predictions of such theories for the system under 
consideration is, therefore, in doubt. 

A more recent theory || extends the dislocation unbinding mechanism for two - 
dimensional melting |22| for systems under external potentials. For a two -dimensional 
triangular solid subjected to an external one -dimensional modulating potential, the only 
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dislocations involved are those which have their Burger's vectors parallel to the troughs of 
the potential. The system, therefore, maps onto an anisotropic, scalar Coulomb gas (or XY 
model) H in contrast to a vector Coulomb gas []22| for the pure 2 — d melting problem. 



Once bound dislocation pairs are integrated out, the melting temperature is obtained as 
a function of the renormalized or "effective" elastic constants which depend on external 
parameters like the strength of the potential, temperature and/or density. Though explicit 
calculations are possible only near the two extreme limits of zero and infinite field intensities 
one can argue effectively that a re -entrant melting transition is expected on general grounds 
quite independent of the detailed nature of the interaction potential for any two -dimensional 
system subject to such external potentials. The actual extent of this region could, of course, 
vary from system to system. In addition, these authors predict that the auto -correlation 
function of the Fourier components of the density (the Debye -Waller correlation function) 
decays algebraically in the solid phase with a universal exponent which depends only on the 
geometry and the magnitude of the reciprocal lattice vector. 

Computer simulation results in this field have so far been inconclusive. Early simula- 
tions []IU[ involving colloidal particles interacting via the Derjaguin, Landau, Verwey and 



Overbeek (DLVO) potential |rjj found a large re -entrant region in apparent agreement with 
later experiments. On closer scrutiny, though, quantitative agreement between simulation 
and experiments on the same system (but with slightly different parameters) appears to be 
poor |J. Subsequent simulations ||iT|-p!3|] have questioned the findings of the earlier com- 
putation and the calculated phase diagram does not show a significant re -entrant liquid 
phase. 

Motivated, in part, by this controversy, we have investigated the freezing/melting be- 
haviour of an unrelated system subjected to similar modulating external potentials. In this 
paper we have computed the phase behaviour of a two dimensional hard disk system in an 
external potential. The pure hard disk system is rather well studied [|i5| -|T8|| by now and 
the nature of the melting transition in the absence of external potentials reasonably well 
explored. Also, there exist colloidal systems with hard interactions [[ljj so that, at least in 



principle, actual experiments using this system are possible. Finally, a hard disk simulation 
is relatively cheap to implement and one can make detailed studies of large systems with- 
out straining computational resources. The main outcome of our calculations, the phase 



iV = 1024 hard disks (diameter a) of density 0.86 < p*(= per 2 ) < 0.91 and the amplitude of 
the external potential < Vq(= /3Vq) < 1000. Within our range of densities, one has a clear 
signature of a re -entrant liquid phase showing that this phenomenon is indeed a general 
one as indicated in Ref. f|. 

The rest of our paper is organized as follows. In Section II we specify the model and the 
simulation method including details of the finite size analysis used. In Section III we present 
our results for the order parameter and its cumulants with a discussion on finite size and 
hysteresis effects. We also present results for the specific heat, order parameter susceptiblity 
and correlation functions which further illustrates the nature of the phase transitions in 
this system. In Section IV we discuss our work in relation to the existing literature on this 
subject, summarize and conclude. 



We study a system of N hard disks of diameter a in a two dimensional box of size S x x Sy 
(S x /Sy = v / 3/2) interacting with the pair potential 4>{rij) between particles % and j with 
distance r^, 



In addition a particle with coordinates (x, y) is exposed to an external periodic potential 



diagram, is shown in Fig. [I]. We have shown results from our simulation of a system of 



II. MODEL AND METHOD 



A. The Model 




(1) 
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of the form: 



V(x iy ) = V Q siYi(2itx/d Q ) (2) 

The constant do in Eq.@ is chosen such that, for a density p = N/S x S y , the modulation is 
commensurate to a triangular lattice of hard disks with nearest neighbor distance a s : do = 
a s \f?>/2 (see Fig. 0). The only parameters which define our system are the reduced density 
per 2 = p* and the reduced potential strength Vq/UbT = V Q *, where ks is the Boltzmann 
constant and T is the temperature. 



B. The Method 



1. Numerical Details 



We perform NVT - Monte Carlo (MC) simulations |T9|j20[| for the system with interactions 
given by Eqs. (|1|) and (|2p for various values of V * and p*. Averages < • > of observables have 
been obtained with the canonical measure. In order to obtain thermodynamic quantities for 
a range of system sizes, we have analyzed various quantities within subsystems as shown in 
Fig. |3|. We have used < • >l to denote averages in subsystems. The subsystems are of size 
L x x L y where L x and L y are chosen as L y = La s and L x = L y \/2>/2 consistent with the 
geometry of the triangular lattice. 

Most of the simulations described below have been done for N = 1024 particles unless 
otherwise indicated. Phase transitions have been studied in most cases by starting in the 
ordered solid and reducing p* for fixed V^*. Runs where the density p* is increased were also 
performed in a few cases. 

A typical simulation run with 4 x 10 7 Monte Carlo steps (MCS) (including 1.5 xlO 7 
MCS for relaxation) took about 50 CPU hours on a PII/500 MHz PC. At high values 
of V * in addition to ordinary (local) MC moves we also used 'through-moves', by which 
particle placements in neighboring troughs are tried. Besides producing faster equilibration, 
including such moves ensures that the formation of dislocations for large V * and p* > 



v^3/2 (do < a) are not artificially hindered since particles can bypass each other — this is 
impossible with purely local MC moves. 

To guarantee good equilibration and averaging, we simulated only systems up to N=1600 
particles in the region of the phase boundary. Systems of N=4096 and N=16384 were 



used only once in Fig. [11], where the interesting region is clearly in the liquid phase and 



equilibration is much easier. 

2. Order parameter 

The nature of the fluid-solid phase transition in two dimensions has been a topic of 



controversy throughout the last forty years |21]-^, [lq , |i^ , ^ , |iS|1 . It is well known that true 
long range positional order is absent in the infinitely large system due to low energy long 
wavelength excitations so that translational correlations decay algebraically. According to 



the dislocation unbinding mechanism |22|J23| the two dimensional solid (with quasi long 



ranged positional order) first melts into a "hexatic" phase with no positional order but 
with quasi long ranged orientational order signified by a non -zero bond orientational order 
parameter ipQ = J^exp(— i68) where 9 is the angle of a bond and the sum is over all distinct 
bonds. A liquid, with no bond orientational order either (ijj G = 0) is produced by a second 
Kosterlitz -Thouless (KT) [^] transition from the hexatic. 

In an external periodic field given by Eq.(Q), however, the bond orientational order 
parameter is non -zero even in the fluid phase |§ . This is because for V * ^ we have now 
a "modulated" liquid, in which local hexagons consisting of the six nearest neighbors of a 
particle are automatically oriented by the external field. Thus < ipe > is non -zero both 
in the (modulated) liquid and the crystalline phase and it cannot be used to study phase 
transitions in this system. The order parameters corresponding to a solid phase are the 
Fourier components of the (non -uniform) density p(r) calculated at the reciprocal lattice 
points {G}. This (infinite) set of numbers are all zero (for G ^ ) in an uniform liquid 
phase and non -zero in a solid. We restrict ourselves to the star consisting of the six smallest 
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reciprocal lattice vectors of the two -dimensional triangular lattice. In the modulated liquid 
phase that is relevant to our system, the Fourier components corresponding to two out of 
these six vectors, viz. those in the direction perpendicular to the troughs of the external 
potential, are non -zero ||. The other four components of this set consisting of those in 
the direction Gi (as defined for the ideal crystal in Fig. |^), and those equivalent to it by 
symmetry, are zero in the (modulated) liquid and non -zero in the solid (if there is true long 
range order). We therefore use the following order parameter: 

N 



^expHG x -rj) 



where fj is the position vector of the j th particle. Note that though the order parameter 
< ipGi > decays to zero with increasing system size in the 2 -d solid — quasi long ranged 
order — this decay, being weak, does not hinder us from distinguishing, in a finite system, 
a modulated liquid from the solid phase with positional order in the G\ direction. 



3. Cumulants 



We have determined phase transition points by the order parameter cumulant intersection 
method 23]. The fourth order cumulant Ul of the order parameter distribution is given by: 



U L (V *,p* 



(3) 



3 < >l 

In case of a continuous transition close to the transition point the cumulant is only a function 
of the ratio of the system size ~ La s and the correlation length £ : Ul(Lcl s /£). Since £ diverges 
at the critical point the cumulants for different system sizes intersect in one point: 17^(0) = 
Ul 2 {0) = U*. Even for first order transitions these cumulants intersect |26| though the value 
U* of Ul at the intersection is not universal any more. The intersection point can, therefore, 
be taken as the phase boundary regardless of the order of the transition. This is useful 
since the order of the melting transition in 2 — d either in the absence pT|- [23|JT6|JT7||2^Jl8 
or with B,[i0|-|i"3|,|9|l external potentials is not unequivocally settled. 



In order to map the phase diagram we systematically vary the system parameters Vq 
and p* to detect order parameter cumulant intersection points which are then identified with 
the phase boundary. It should be noted that though the order parameter (defined for long 



range positional order) vanishes f22| , |23|| with increasing system size in the crystalline phase, 
its cumulants are well defined and can be used to determine phase boundaries. For large 
L the cumulants approach the value 2/3 in the solid phase and 1/3 in the liquid [|30| so 
that they are guaranteed to intersect ! For very large V^j* we do not find an unique point 
of intersection for Ul, instead the cumulants for various values of L collapse onto a single 
curve. In this case the onset of the collapse is taken as the "intersection" density. It is 



curious to note that this behaviour is, in fact, typical of the anisotropic XY model |]27 



In this case although the order parameter cumulants have an intersection point, the value 
of the cumulant at the intersection differs for various anisotropies and drifts towards a 
limiting value at zero anisotropy. The intersection "point" therefore changes to a "line" of 
intersections for different system sizes and for small anisotropies. In our system, for the 
large Vq we see similar behaviour. 



III. RESULTS AND DISCUSSION 

A. Order parameter and cumulants 

In Fig. f| we present data for the average order parameter < ipQ x > and its cumulants as 
functions of the density for V^* = 0.05 and V^* = 0.5 calculated within various subsystems. 
In both cases < ipGi >l and Ul increase with p* with a sharpening of the structures for 
increasing L. As discussed above, we observe that for any density increasing subsystem size 
L depresses < >l- The cumulant Ul, on the other hand, approaches limiting values 
(2/3 for solid and 1/3 for liquid). The values of the cumulants are higher for larger L in 
the ordered (solid) phase and vice versa in the disordered (modulated liquid) phase thus 
resulting in an intersection point — the transition density. 
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In Fig. [| we compare the density dependence of the average order parameter < if)a 1 > 
(calculated over the entire system) for different V * and for the same system size (N = 1024). 
With increasing V * the turning point in < ipQ 1 > (p*) is shifted to lower densities and then 
for even larger V * values to higher densities. This indicates, already, that the system prefers 
having smaller transition densities for intermediate values of Vq compared to smaller and 
higher V * values — i.e. we have a re -entrant transition. 

In Fig. H we show a systematic study of < i/;g 1 >l and {/l as a function of V * at the 
density p* = 0.89 for different L- values. Maxima in the < ipGi >l- and Ul- curves are 
found near V * = 2. Again we note that the < tpd >ir values decrease with increasing L 
(see Fig. §(a)). The cumulants Ul, on the other hand increases with L for intermediate 
values of V * (the ordered, solid phase) and decreases with L for either large or small V * (the 
disordered, liquid phase) resulting in intersection points indicating two consecutive phase 
transitions (see Fig. 0(b)). 

If V * is increased to 20 the value of the cumulant at intersection U* is shifted upwards, 
see Fig. 0(a). For very high Vq*- values the cumulant curves for different L merge on the 
high density side, see Fig. 0(b) (see discussion in Section [11 B 3[ ). In Fig. || the cumulant 



intersection values are shown as a function of V *, where for large VQ*-values the value at 
the onset of the merging is shown. We observe that U* is not an universal number but, 
nevertheless, goes to a limiting value for large V * [27| . 



In Fig. | we show < ip Gl > as a function of the density with V * = 0.5 for different 
iV-values. The general features of < ijjoi > as discussed above is retained though there is a 
shift of the turning point to slightly higher densities with increasing N. The effect on the 
phase diagram is discussed in Section |III Q . 

B. Susceptibility, specific heat, finite size effects and hysteresis 

In addition to < ipG 1 >l and Ul we have computed the order parameter susceptibility 
Xgi and the specific heat for different system- and subsystem-sizes. 
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The order parameter susceptibility XGi is defined as |3(J: 



k B Tx Gl = L 2 [<(^ Gl ) 2 > - & Gl f] (4) 

In Fig. 0(a) we show xg ± as a function of p* at Vq* = 0.05 for different L- values. The 
increase of Xd with increasing L signals the presence of a phase transition in the density 
range where the transition has been found by cumulant intersection techniques {pi m 0.896). 
In Fig. [l^(b) XGi is shown for the same system size (N = 1024) and various V^*-values. We 
note that the density of the XGi~ maxima are smallest for the intermediate value of Vq 
which again show that for these V^*- values the transition density is lowest. Compared to 
the cumulant intersection values, Xd maxima are located at slightly smaller densities (see 
also Section [III C| ) which may be due to finite size effects, which often show the feature that 
phase transition points in finite systems are shifted to slightly different values depending on 
the observable under investigation. In particular one expects (and we get) a shift towards 
parameter values in the disordered region (here a liquid, i.e. low densities) for the order 
parameter and the susceptibility as compared to the cumulant intersection parameters. 

We have also calculated the specific heat Cy (Fig. [TTj) as a function of the density for 
Vq = 0.2 with N = 4096 and iV = 16384. For a second order transition, the maximum of the 
specific heat scales with the system size as Cy ax ~ L a ' v where a and v are critical exponents. 
For a first order transition, on the other hand Cy ax ~ L d where d is the dimensionality (= 2 
in our case). We, however, do not see any of this behaviour. In contrast, the specific heat 
is relatively featureless. Although it shows a peak, surprisingly, the height of this peak is 
almost insensitive to system size. This is a strong indication that the phase transition we 
observe is unconventional and is KT like |22||| . Further, as expected for such transitions, the 
maximum does not lie at the density where the cumulants intersect and it would be incorrect 
to identify specific heat maxima with the phase boundary (see discussion in Section |TV| ). 

In order to study the effect of the path taken through the parameter space on the 
location of the cumulant intersection densities, we compared Ul as a function of the density 
for Vq = 0.5 as obtained from two runs. In Fig. f|(d) we have already presented the data 
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for a run where the density was decreased systematically. In Fig. [12| we present data for 
a second run where the density was increased instead. We find negligibly small hysteresis 
effects on the cumulants as well as on the value of intersection density. This shows that the 
transition points are not affected much by the path through the parameter space. 



C. The Phase Diagram 

We have obtained the phase diagram of the system for .86 < p* < .91 and < Vq < 1000. 
For each density and Vq - value we computed cumulants Ul for a range of subsystem sizes 
L and located intersection points which we identify with the phase boundary. The resulting 
locus of the transition points is shown in the phase diagram, see Fig. [[]. At very small 
values of V^* we find good agreement of our transition densities with the melting densities 
(p m pa 0.91) known from literature |17| on the pure hard disk solid (Vq = 0). The values 
of the transition density initially drop and subsequently rise as Vq increases. The minimum 
transition density is found for V^* ~ 1 — 2. These transition points separate a high density 
solid from a low density modulated liquid. Thus, at a properly chosen density, we observe 
an initial freezing transition followed by a re -entrant melting at a higher V^*- value. Such 
an effect had been found earlier in experiments on colloidal systems in an external laser 
field [If. 

In order to quantify finite size effects on the phase diagram, we have computed the 
transition points for different system sizes. The resulting phase diagrams are shown in 



Fig. [L3|. We note that due to residual finite size effects with increasing system size all 
transition points are slightly shifted to higher densities, the structure of the phase diagram 
with a pronounced minimum at intermediate values of V^* is not affected by this shift. 



D. Correlation functions 

The Debye- Waller- correlation function is defined as follows: 
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C di (R) =< e *GiWR)-m) > 

where R points to the elementary cell of the ideal lattice, and u{R) is the deviation of the 
actual particle position from the ideal lattice: r = R + u(R) . In this case we have chosen 
the direction of R to lie along the y axis (i.e. along the troughs of the potential). 

We have also computed the spatial correlation function g(y) which is the pair correlation 
function in the y-direction. We compute it in the following way: for a particle i, g(y) dy oc 
number of particles j for which: \y { - yj\ G [y, y + dy] and | Ob ^ ^ j | 

do/2, normalized so 

that g(y) — > 1 as y — > oo. 



These correlation functions are plotted in Fig. [TJ] as functions of y. The Debye- Waller 
correlation function C(j 1 (y) and the correlation function g(y) along the potential valley 
are compared in Fig. |l4|(a) at a density just below the transition. We see that the decay 
of the maxima of g(y) as function of y is similar to the decay in Ca 1 (y). The decay of 
C*Gi(y) is analyzed in more detail in Fig. 0(b) for parameter values in the liquid and in 
the solid phase. In the liquid phase the decay is exponential while in the solid region it is 
algebraic: C^iy) ~ y~ m ±. Taking the data-points in the crystal which are closest to the 
phase boundary for each V *, we get r] Gl in the range of 0.20 . . . 0.27. The exponent i] Gl is 
predicted to be universal and equal to 1/4 for our geometry, so this value is consistent 
with our numerical results. 

E. Scaling behavior 

We next try to determine the order of the phase transitions encountered in this system 
for various values of V^*. In order to investigate this issue we studied the scaling behavior 
of the order parameter, susceptibility and the order parameter cumulant near the phase 
boundary for a small (0.5) and a large (1000) Vq. From finite size scaling theory (for an 



overview see Ref. [p0| ) we expect these quantities to scale as [[28[ : 
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(5) 



XL k B TL- c ~ p(L/0 (6) 

t4 ~ WO (7) 

Here b = /3/v, c = ^/v (for critical scaling) and /, g>, /i are scaling functions. The correlation 
length £ diverges as £ oc (1 — p/ p c )~ v for an ordinary critical point, while for a KT- transition 
we have an essential singularity and £ oc exp(a(l — p/ p c )~ v ). 

According to general arguments given in Ref. 0, we expect that for a finite lattice, the 
identification of the properties of our system with those of the anisotropic XY model should 
improve with increasing V *. Indeed, for large V *, scaling according to the KT- theory seems 
to be supported by our data. In Fig. (|T5| ) we have plotted the left hand sides of Eqs.(^) 
and (H) versus L/£ for V * = 1000, where data points for 0.86 < p* < 0.898 have been 
considered and p* = 0.902, obtained by cumulant intersection. In order not to introduce 
an unwarranted bias, we have separately considered (a) ordinary critical scaling and (b) 
a KT scaling forms and adjusted the values of the parameters b, c and v till we obtained 
collapse of our data onto a single curve determined by a least square estimator. Though 
good collapse of our data is observed both in (a) and (b), the numerical values for u, rj = lb 
and c = 2 — 7] for KT scaling (b « 0.138, c ~ 1.70, v w 0.44) are close to the predicted 
values = r]/2 = 1/8, c = 1.75, z> = 0.5). 

The situation is less straight -forward for Vq = 0.5. The critical parameters were obtained 



in this case for densities 0.85 < p* < 0.876, with p* = 0.878. In Fig. [16] (a) the data collapse 
looks slightly better than in (b), such that relying on this data alone one may conclude 
that KT- scaling in this region of the phase diagram seems less likely. It must be kept in 
mind though that for small values of V^* in a finite system, the analysis of the data would 
be complicated by crossover effects. Strictly for Vq = we do not have a correspondence 
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with the anisotropic XY model but rather with a vector Coulomb gas with a different 
set of exponents. Our results for the numerical values of the parameters are summarized in 
Table I. 

In summary, from the scaling analysis in Figs. ([16]) and ( [T5|) a KT- scenario at least for 
large V^*- values seems likely. This is supported by the behavior of the cumulants as well 



(see Sect. [Ill Aj) . A more precise classification of the phase transitions with the present data 
and system sizes is not easy. This topic is left for future work, in particular we plan to 
compute the elastic properties of the system by a method recently developed for the hard- 
disk system 0j2^] and to test the KT- predictions 0. 



IV. SUMMARY AND CONCLUSION 

In summary, we have calculated the phase diagram of a two dimensional system of hard 
disks in an external sinusoidal potential. We found freezing followed by re -entrant melting 
transitions over a significant region of the phase diagram in tune with previous experiments 
on colloids @-§| and with the expectations of a recent dislocation unbinding theory ||. One 
of the main features of our calculation is the method used to locate phase boundaries. In 
contrast to earlier simulations []T0|-[T3|] which used either the jump of the order parameter, or 
specific heat maxima to locate the phase transition, we have used the more reliable cumulant 
intersection method. It must be noted that the specific heat in this system does not show 
a strong peak at the phase transition density so that its use may lead to confusing results. 
This, in our opinion, may be the reason for part of the controversy in this field. It is possible 
that earlier simulations which used smaller systems and no systematic finite size analysis may 
have overlooked this feature of Cy which becomes apparent only in computations involving 
large system sizes. We have shown that finite size scaling of the order parameter cumulants 
as obtained from sub-system or sub-block analysis, on the other hand, yields an accurate 
phase diagram. 

What is the order of the phase transition? We know that JTBHTBl for the pure hard 
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disk system in two -dimensions this question is quite difficult to answer and our present 



understanding jlB) is that hard disks in two dimensions show a KTHNY transition. This 
transition, however, lies very close (in an appropriately expanded parameter space) to a first 
order boundary so that crossover effects may be significant. The present system has an im- 
posed external periodic potential which stabilizes the hexatic phase ]9[] and an (anisotropic) 
KT transition |J is expected. Our results show several features which suggest that this is, 
perhaps, what we have. Though we have discussed these observations in the rest of the 
paper, we list them below for clarity: 

• The behaviour of the value of the cumulants at intersection U* is similar to an earlier 



work |27| on the anisotropic XY system which shows a KT transition. 



• The specific heat is relatively featureless and does not scale with system size in a 
fashion expected of a true first order or conventional continuous transition. 

• The decay of the correlation functions is similar to what is predicted || for an 
anisotropic scalar Coulomb gas. 

• For large Vq- values the scaling of the order parameter, the susceptibility and the 
cumulant may be described by the KT- theory. 

Of course, in order to resolve this issue unambiguously yet larger simulations are required. 
Also, we need to compute elastic properties |2^Jl8|| of this system in order to compare directly 



with the results of Ref. ||. Work along these lines is in progress. 

Before we end we would like to point out that after completion of this work and prior 
to the submission of this manuscript we received a preprint where the same system as 
ours has been studied using simulations. The phase diagram obtained by these authors is 
similar to ours (thus confirming and corroborating our results), though there exist some 
quantitative differences. These differences may be attributed to the absence of systematic 
finite size scaling in the latter work. 
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lowing probability distribution in the (modulated) liquid and L ^> £ [pG 



k B Tx 



k B Tx 



which yields: Ul = 1/3 , < tpQ >= fcj ^l* , and x = x(l _ 7r /4) for x defined as in (f|). 
Because ^ is very similar in spirit, the same results apply for w(ipe) in a isotropic liquid 
and L ^> £. 



19 



V* 


b 


c 


V 


b 


c 


V 


a 


1000 


0.13(1) 


1.68(5) 


2.25(25) 


0.138(8) 


1.70(5) 


0.44(3) 


1.05(25) 


0.5 


0.152(5) 


1.65(6) 


1.06(13) 


0.170(12) 


1.83(4) 


0.38(10) 


1.0(2) 



Table I Parameters in the scaling plots (Figs. dl6|) and (|15D) for V * = 0.5 and Vq = 1000. 
The first three parameter columns are for critical scaling, the last four for KT- scaling. 
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FIGURES 




100 1000 



FIG. 1. Phase diagram in the p* / Vq- plane. Transition points for transitions from the solid 
to the modulated liquid have been obtained by the order parameter cumulant intersection method. 
In order to map the phase diagram we scanned in p* for every Vq, starting from the high density 
(solid) region. The system size is N = 1024. 
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FIG. 2. Schematic picture of the system geometry showing the direction G\ along which 
crystalline order develops in the modulated liquid. The four vectors obtained by rotating G\ 
anti -clockwise by 60° and/or reflecting about the origin are equivalent. The parameters cLq and a s 
are also shown. The size of the box is S x x S y and the modulating potential is V . 
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FIG. 3. Schematic picture showing subboxes of size L (here L = 3) used in the finite size 
scaling analysis (see text). 
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FIG. 4. Order parameter < tpQ 1 >l ((a) and (c)) and order parameter cumulant Ul ((b) and 
(d)) versus p* in subsystems of size L for reduced potential amplitudes Vg" = 0.05 ((a) and (b)) 
and Vq = 0.5 ((c) and (d)), N = 1024. Unless otherwise stated, lines connecting data points in 
this and the rest of the figures are for visual guidance. 
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FIG. 6. Order parameter < tpa 1 >l (a) and order parameter cumulant Ul (b) versus Vq at 
a constant density p* = 0.89 for different L, the system size is N = 1024. 
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FIG. 8. Values U* of the order parameter cumulants at the intersection points versus V *. 
The shown data at the largest four values of V * are taken at the onset of the cumulants curves 
merging (see text). The system size is N = 1024. 
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p p 

FIG. 10. Order parameter susceptibilities versus density for: (a) constant Vq = 0.05 and 
different L values, N = 1024, (b) full system size (N = 1024) and different Vq - values. 
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FIG. 11. Specific heat per particle versus density at constant Vq = 0.2 and different system 
sizes (TV = 4096, 16384). 
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0.89 



FIG. 12. Order parameter cumulant Ul versus density at constant Vq = 0.5 and different L, 
N = 1024. Values are obtained by successively compressing the system from one density to the 
next higher density. For a corresponding picture obtained by successively expanding the system 
see Fig. |(d). 
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FIG. 13. Phase diagram in the p* / V *- plane for different system sizes (N = 400, 676, 1024, 
1600). 
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FIG. 14. Debye- Waller (Cdiy)) correlation- and pair distribution- (g(y)) functions as func- 
tions of y parallel to the potential minima for fixed p* = 0.86 and Vq = 2 (a) and Debye- Waller 
correlation function versus y for fixed p* = 0.88 and different Vq = 0.1,2, 1000 (b). Lines in the 
upper right corner of (b) show the functions f(y) = y~ w (w = 0.1, 0.2, 0.25) for comparison. The 
system size is N = 1024. 
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FIG. 15. Scaling plots for the order parameter cumulant and the susceptibility (inset) for 
Vq = 1000 assuming (a) critical and (b) KT- scaling. The system size is N = 1024, for £ we have 
used the expressions given after Eq.(Q). 
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FIG. 16. Scaling plots for the order parameter and the susceptibility (inset) for Vq = 0.5 
assuming (a) critical and (b) KT- scaling. The system size is N = 1024, for £ we have used the 
expressions given after Eq.(^). 
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